notebook.community

Edit and run 





In [1]:



    


import numpy as np
%matplotlib inline
import matplotlib.pyplot as plt

from sympy import *
from mpmath import degrees, radians

init_printing()

This derivation is based on the USGS Professional Paper 1395: Map Projections — A Working Manual, ch 14, pp 98 – 103.





In [2]:



    


radius = Symbol('R')
semimajor, ellipticity = symbols('a e', positive=True)
standard_parallel_1, standard_parallel_2 = symbols('phi_1 phi_2')
central_lat, central_lon = symbols('phi_0 lambda_0')
point_lat, point_lon = symbols('phi lambda')

radius, semimajor, ellipticity, \
standard_parallel_1, standard_parallel_2, \
central_lat, central_lon, \
point_lat, point_lon



    


















    
Out[2]:







$$\left ( R, \quad a, \quad e, \quad \phi_{1}, \quad \phi_{2}, \quad \phi_{0}, \quad \lambda_{0}, \quad \phi, \quad \lambda\right )$$
















In [3]:



    


def subdict(eq, vals):
    for k, v in vals.items():
        eq = eq.subs(k, v)
    return eq

Formulas for the Sphere

NOTE: Formulas are defined in reverse (mostly) from the paper, so that they can build on each other.





In [4]:



    


# Eq. 14-6
ns = (sin(standard_parallel_1) + sin(standard_parallel_2)) / 2
ns



    


















    
Out[4]:







$$\frac{1}{2} \sin{\left (\phi_{1} \right )} + \frac{1}{2} \sin{\left (\phi_{2} \right )}$$
















In [5]:



    


# Eq. 14-5
Cs = cos(standard_parallel_1)**2 + 2 * ns * sin(standard_parallel_1)
Cs



    


















    
Out[5]:







$$\left(\sin{\left (\phi_{1} \right )} + \sin{\left (\phi_{2} \right )}\right) \sin{\left (\phi_{1} \right )} + \cos^{2}{\left (\phi_{1} \right )}$$
















In [6]:



    


# Eq. 14-3a
rho0s = radius * (Cs - 2 * ns * sin(central_lat))**0.5 / ns
rho0s



    


















    
Out[6]:







$$\frac{R}{\frac{1}{2} \sin{\left (\phi_{1} \right )} + \frac{1}{2} \sin{\left (\phi_{2} \right )}} \left(- \left(\sin{\left (\phi_{1} \right )} + \sin{\left (\phi_{2} \right )}\right) \sin{\left (\phi_{0} \right )} + \left(\sin{\left (\phi_{1} \right )} + \sin{\left (\phi_{2} \right )}\right) \sin{\left (\phi_{1} \right )} + \cos^{2}{\left (\phi_{1} \right )}\right)^{0.5}$$
















In [7]:



    


# Eq. 14-4
thetas = ns * (point_lon - central_lon)
thetas



    


















    
Out[7]:







$$\left(\lambda - \lambda_{0}\right) \left(\frac{1}{2} \sin{\left (\phi_{1} \right )} + \frac{1}{2} \sin{\left (\phi_{2} \right )}\right)$$
















In [8]:



    


# Eq. 14-3
rhos = radius * (Cs - 2 * ns * sin(point_lat))**0.5 / ns
rhos



    


















    
Out[8]:







$$\frac{R}{\frac{1}{2} \sin{\left (\phi_{1} \right )} + \frac{1}{2} \sin{\left (\phi_{2} \right )}} \left(- \left(\sin{\left (\phi_{1} \right )} + \sin{\left (\phi_{2} \right )}\right) \sin{\left (\phi \right )} + \left(\sin{\left (\phi_{1} \right )} + \sin{\left (\phi_{2} \right )}\right) \sin{\left (\phi_{1} \right )} + \cos^{2}{\left (\phi_{1} \right )}\right)^{0.5}$$
















In [9]:



    


# Eq. 14-1
xs = rhos * sin(thetas)
# Eq. 14-2
ys = rho0s - rhos * cos(thetas)



    











In [10]:



    


xs



    


















    
Out[10]:







$$\frac{R \sin{\left (\left(\lambda - \lambda_{0}\right) \left(\frac{1}{2} \sin{\left (\phi_{1} \right )} + \frac{1}{2} \sin{\left (\phi_{2} \right )}\right) \right )}}{\frac{1}{2} \sin{\left (\phi_{1} \right )} + \frac{1}{2} \sin{\left (\phi_{2} \right )}} \left(- \left(\sin{\left (\phi_{1} \right )} + \sin{\left (\phi_{2} \right )}\right) \sin{\left (\phi \right )} + \left(\sin{\left (\phi_{1} \right )} + \sin{\left (\phi_{2} \right )}\right) \sin{\left (\phi_{1} \right )} + \cos^{2}{\left (\phi_{1} \right )}\right)^{0.5}$$
















In [11]:



    


ys



    


















    
Out[11]:







$$- \frac{R \cos{\left (\left(\lambda - \lambda_{0}\right) \left(\frac{1}{2} \sin{\left (\phi_{1} \right )} + \frac{1}{2} \sin{\left (\phi_{2} \right )}\right) \right )}}{\frac{1}{2} \sin{\left (\phi_{1} \right )} + \frac{1}{2} \sin{\left (\phi_{2} \right )}} \left(- \left(\sin{\left (\phi_{1} \right )} + \sin{\left (\phi_{2} \right )}\right) \sin{\left (\phi \right )} + \left(\sin{\left (\phi_{1} \right )} + \sin{\left (\phi_{2} \right )}\right) \sin{\left (\phi_{1} \right )} + \cos^{2}{\left (\phi_{1} \right )}\right)^{0.5} + \frac{R}{\frac{1}{2} \sin{\left (\phi_{1} \right )} + \frac{1}{2} \sin{\left (\phi_{2} \right )}} \left(- \left(\sin{\left (\phi_{1} \right )} + \sin{\left (\phi_{2} \right )}\right) \sin{\left (\phi_{0} \right )} + \left(\sin{\left (\phi_{1} \right )} + \sin{\left (\phi_{2} \right )}\right) \sin{\left (\phi_{1} \right )} + \cos^{2}{\left (\phi_{1} \right )}\right)^{0.5}$$

Formulas for the Ellipsoid





In [12]:



    


# Eq. 14-15
m = cos(point_lat) / (1 - ellipticity**2 * sin(point_lat)**2)**0.5
m1 = cos(standard_parallel_1) / (1 - ellipticity**2 * sin(standard_parallel_1)**2)**0.5
m2 = cos(standard_parallel_2) / (1 - ellipticity**2 * sin(standard_parallel_2)**2)**0.5
m, m1, m2



    


















    
Out[12]:







$$\left ( \frac{\cos{\left (\phi \right )}}{\left(- e^{2} \sin^{2}{\left (\phi \right )} + 1\right)^{0.5}}, \quad \frac{\cos{\left (\phi_{1} \right )}}{\left(- e^{2} \sin^{2}{\left (\phi_{1} \right )} + 1\right)^{0.5}}, \quad \frac{\cos{\left (\phi_{2} \right )}}{\left(- e^{2} \sin^{2}{\left (\phi_{2} \right )} + 1\right)^{0.5}}\right )$$
















In [13]:



    


# Eq 3-12
q = (1 - ellipticity**2) * (
    sin(point_lat) / (1 - ellipticity**2 * sin(point_lat)**2) -
    1 / (2 * ellipticity) * log(
        (1 - ellipticity * sin(point_lat)) /
        (1 + ellipticity * sin(point_lat))
        )
    )
q0 = (1 - ellipticity**2) * (
    sin(central_lat) / (1 - ellipticity**2 * sin(central_lat)**2) -
    1 / (2 * ellipticity) * log(
        (1 - ellipticity * sin(central_lat)) /
        (1 + ellipticity * sin(central_lat))
        )
    )
q1 = (1 - ellipticity**2) * (
    sin(standard_parallel_1) / (1 - ellipticity**2 * sin(standard_parallel_1)**2) -
    1 / (2 * ellipticity) * log(
        (1 - ellipticity * sin(standard_parallel_1)) /
        (1 + ellipticity * sin(standard_parallel_1))
        )
    )
q2 = (1 - ellipticity**2) * (
    sin(standard_parallel_2) / (1 - ellipticity**2 * sin(standard_parallel_2)**2) -
    1 / (2 * ellipticity) * log(
        (1 - ellipticity * sin(standard_parallel_2)) /
        (1 + ellipticity * sin(standard_parallel_2))
        )
    )

q, q0, q1, q2



    


















    
Out[13]:







$$\left ( \left(- e^{2} + 1\right) \left(\frac{\sin{\left (\phi \right )}}{- e^{2} \sin^{2}{\left (\phi \right )} + 1} - \frac{1}{2 e} \log{\left (\frac{- e \sin{\left (\phi \right )} + 1}{e \sin{\left (\phi \right )} + 1} \right )}\right), \quad \left(- e^{2} + 1\right) \left(\frac{\sin{\left (\phi_{0} \right )}}{- e^{2} \sin^{2}{\left (\phi_{0} \right )} + 1} - \frac{1}{2 e} \log{\left (\frac{- e \sin{\left (\phi_{0} \right )} + 1}{e \sin{\left (\phi_{0} \right )} + 1} \right )}\right), \quad \left(- e^{2} + 1\right) \left(\frac{\sin{\left (\phi_{1} \right )}}{- e^{2} \sin^{2}{\left (\phi_{1} \right )} + 1} - \frac{1}{2 e} \log{\left (\frac{- e \sin{\left (\phi_{1} \right )} + 1}{e \sin{\left (\phi_{1} \right )} + 1} \right )}\right), \quad \left(- e^{2} + 1\right) \left(\frac{\sin{\left (\phi_{2} \right )}}{- e^{2} \sin^{2}{\left (\phi_{2} \right )} + 1} - \frac{1}{2 e} \log{\left (\frac{- e \sin{\left (\phi_{2} \right )} + 1}{e \sin{\left (\phi_{2} \right )} + 1} \right )}\right)\right )$$
















In [14]:



    


# Eq. 14-14
ne = (m1**2 - m2**2) / (q2 - q1)
ne



    


















    
Out[14]:







$$\frac{\frac{\cos^{2}{\left (\phi_{1} \right )}}{\left(- e^{2} \sin^{2}{\left (\phi_{1} \right )} + 1\right)^{1.0}} - \frac{\cos^{2}{\left (\phi_{2} \right )}}{\left(- e^{2} \sin^{2}{\left (\phi_{2} \right )} + 1\right)^{1.0}}}{- \left(- e^{2} + 1\right) \left(\frac{\sin{\left (\phi_{1} \right )}}{- e^{2} \sin^{2}{\left (\phi_{1} \right )} + 1} - \frac{1}{2 e} \log{\left (\frac{- e \sin{\left (\phi_{1} \right )} + 1}{e \sin{\left (\phi_{1} \right )} + 1} \right )}\right) + \left(- e^{2} + 1\right) \left(\frac{\sin{\left (\phi_{2} \right )}}{- e^{2} \sin^{2}{\left (\phi_{2} \right )} + 1} - \frac{1}{2 e} \log{\left (\frac{- e \sin{\left (\phi_{2} \right )} + 1}{e \sin{\left (\phi_{2} \right )} + 1} \right )}\right)}$$
















In [15]:



    


# Eq. 14-13
Ce = m1**2 + ne * q1
Ce



    


















    
Out[15]:







$$\frac{\left(- e^{2} + 1\right) \left(\frac{\sin{\left (\phi_{1} \right )}}{- e^{2} \sin^{2}{\left (\phi_{1} \right )} + 1} - \frac{1}{2 e} \log{\left (\frac{- e \sin{\left (\phi_{1} \right )} + 1}{e \sin{\left (\phi_{1} \right )} + 1} \right )}\right) \left(\frac{\cos^{2}{\left (\phi_{1} \right )}}{\left(- e^{2} \sin^{2}{\left (\phi_{1} \right )} + 1\right)^{1.0}} - \frac{\cos^{2}{\left (\phi_{2} \right )}}{\left(- e^{2} \sin^{2}{\left (\phi_{2} \right )} + 1\right)^{1.0}}\right)}{- \left(- e^{2} + 1\right) \left(\frac{\sin{\left (\phi_{1} \right )}}{- e^{2} \sin^{2}{\left (\phi_{1} \right )} + 1} - \frac{1}{2 e} \log{\left (\frac{- e \sin{\left (\phi_{1} \right )} + 1}{e \sin{\left (\phi_{1} \right )} + 1} \right )}\right) + \left(- e^{2} + 1\right) \left(\frac{\sin{\left (\phi_{2} \right )}}{- e^{2} \sin^{2}{\left (\phi_{2} \right )} + 1} - \frac{1}{2 e} \log{\left (\frac{- e \sin{\left (\phi_{2} \right )} + 1}{e \sin{\left (\phi_{2} \right )} + 1} \right )}\right)} + \frac{\cos^{2}{\left (\phi_{1} \right )}}{\left(- e^{2} \sin^{2}{\left (\phi_{1} \right )} + 1\right)^{1.0}}$$
















In [16]:



    


# Eq. 14-12a
rho0e = semimajor * (Ce - ne * q0)**0.5 / ne
rho0e



    


















    
Out[16]:







$$\frac{a}{\frac{\cos^{2}{\left (\phi_{1} \right )}}{\left(- e^{2} \sin^{2}{\left (\phi_{1} \right )} + 1\right)^{1.0}} - \frac{\cos^{2}{\left (\phi_{2} \right )}}{\left(- e^{2} \sin^{2}{\left (\phi_{2} \right )} + 1\right)^{1.0}}} \left(- \left(- e^{2} + 1\right) \left(\frac{\sin{\left (\phi_{1} \right )}}{- e^{2} \sin^{2}{\left (\phi_{1} \right )} + 1} - \frac{1}{2 e} \log{\left (\frac{- e \sin{\left (\phi_{1} \right )} + 1}{e \sin{\left (\phi_{1} \right )} + 1} \right )}\right) + \left(- e^{2} + 1\right) \left(\frac{\sin{\left (\phi_{2} \right )}}{- e^{2} \sin^{2}{\left (\phi_{2} \right )} + 1} - \frac{1}{2 e} \log{\left (\frac{- e \sin{\left (\phi_{2} \right )} + 1}{e \sin{\left (\phi_{2} \right )} + 1} \right )}\right)\right) \left(- \frac{\left(- e^{2} + 1\right) \left(\frac{\sin{\left (\phi_{0} \right )}}{- e^{2} \sin^{2}{\left (\phi_{0} \right )} + 1} - \frac{1}{2 e} \log{\left (\frac{- e \sin{\left (\phi_{0} \right )} + 1}{e \sin{\left (\phi_{0} \right )} + 1} \right )}\right) \left(\frac{\cos^{2}{\left (\phi_{1} \right )}}{\left(- e^{2} \sin^{2}{\left (\phi_{1} \right )} + 1\right)^{1.0}} - \frac{\cos^{2}{\left (\phi_{2} \right )}}{\left(- e^{2} \sin^{2}{\left (\phi_{2} \right )} + 1\right)^{1.0}}\right)}{- \left(- e^{2} + 1\right) \left(\frac{\sin{\left (\phi_{1} \right )}}{- e^{2} \sin^{2}{\left (\phi_{1} \right )} + 1} - \frac{1}{2 e} \log{\left (\frac{- e \sin{\left (\phi_{1} \right )} + 1}{e \sin{\left (\phi_{1} \right )} + 1} \right )}\right) + \left(- e^{2} + 1\right) \left(\frac{\sin{\left (\phi_{2} \right )}}{- e^{2} \sin^{2}{\left (\phi_{2} \right )} + 1} - \frac{1}{2 e} \log{\left (\frac{- e \sin{\left (\phi_{2} \right )} + 1}{e \sin{\left (\phi_{2} \right )} + 1} \right )}\right)} + \frac{\left(- e^{2} + 1\right) \left(\frac{\sin{\left (\phi_{1} \right )}}{- e^{2} \sin^{2}{\left (\phi_{1} \right )} + 1} - \frac{1}{2 e} \log{\left (\frac{- e \sin{\left (\phi_{1} \right )} + 1}{e \sin{\left (\phi_{1} \right )} + 1} \right )}\right) \left(\frac{\cos^{2}{\left (\phi_{1} \right )}}{\left(- e^{2} \sin^{2}{\left (\phi_{1} \right )} + 1\right)^{1.0}} - \frac{\cos^{2}{\left (\phi_{2} \right )}}{\left(- e^{2} \sin^{2}{\left (\phi_{2} \right )} + 1\right)^{1.0}}\right)}{- \left(- e^{2} + 1\right) \left(\frac{\sin{\left (\phi_{1} \right )}}{- e^{2} \sin^{2}{\left (\phi_{1} \right )} + 1} - \frac{1}{2 e} \log{\left (\frac{- e \sin{\left (\phi_{1} \right )} + 1}{e \sin{\left (\phi_{1} \right )} + 1} \right )}\right) + \left(- e^{2} + 1\right) \left(\frac{\sin{\left (\phi_{2} \right )}}{- e^{2} \sin^{2}{\left (\phi_{2} \right )} + 1} - \frac{1}{2 e} \log{\left (\frac{- e \sin{\left (\phi_{2} \right )} + 1}{e \sin{\left (\phi_{2} \right )} + 1} \right )}\right)} + \frac{\cos^{2}{\left (\phi_{1} \right )}}{\left(- e^{2} \sin^{2}{\left (\phi_{1} \right )} + 1\right)^{1.0}}\right)^{0.5}$$
















In [17]:



    


# Eq. 14-4
thetae = ne * (point_lon - central_lon)
thetae



    


















    
Out[17]:







$$\frac{\left(\lambda - \lambda_{0}\right) \left(\frac{\cos^{2}{\left (\phi_{1} \right )}}{\left(- e^{2} \sin^{2}{\left (\phi_{1} \right )} + 1\right)^{1.0}} - \frac{\cos^{2}{\left (\phi_{2} \right )}}{\left(- e^{2} \sin^{2}{\left (\phi_{2} \right )} + 1\right)^{1.0}}\right)}{- \left(- e^{2} + 1\right) \left(\frac{\sin{\left (\phi_{1} \right )}}{- e^{2} \sin^{2}{\left (\phi_{1} \right )} + 1} - \frac{1}{2 e} \log{\left (\frac{- e \sin{\left (\phi_{1} \right )} + 1}{e \sin{\left (\phi_{1} \right )} + 1} \right )}\right) + \left(- e^{2} + 1\right) \left(\frac{\sin{\left (\phi_{2} \right )}}{- e^{2} \sin^{2}{\left (\phi_{2} \right )} + 1} - \frac{1}{2 e} \log{\left (\frac{- e \sin{\left (\phi_{2} \right )} + 1}{e \sin{\left (\phi_{2} \right )} + 1} \right )}\right)}$$
















In [18]:



    


# Eq. 14-12
rhoe = semimajor * (Ce - ne * q)**0.5 / ne
rhoe



    


















    
Out[18]:







$$\frac{a}{\frac{\cos^{2}{\left (\phi_{1} \right )}}{\left(- e^{2} \sin^{2}{\left (\phi_{1} \right )} + 1\right)^{1.0}} - \frac{\cos^{2}{\left (\phi_{2} \right )}}{\left(- e^{2} \sin^{2}{\left (\phi_{2} \right )} + 1\right)^{1.0}}} \left(- \left(- e^{2} + 1\right) \left(\frac{\sin{\left (\phi_{1} \right )}}{- e^{2} \sin^{2}{\left (\phi_{1} \right )} + 1} - \frac{1}{2 e} \log{\left (\frac{- e \sin{\left (\phi_{1} \right )} + 1}{e \sin{\left (\phi_{1} \right )} + 1} \right )}\right) + \left(- e^{2} + 1\right) \left(\frac{\sin{\left (\phi_{2} \right )}}{- e^{2} \sin^{2}{\left (\phi_{2} \right )} + 1} - \frac{1}{2 e} \log{\left (\frac{- e \sin{\left (\phi_{2} \right )} + 1}{e \sin{\left (\phi_{2} \right )} + 1} \right )}\right)\right) \left(- \frac{\left(- e^{2} + 1\right) \left(\frac{\sin{\left (\phi \right )}}{- e^{2} \sin^{2}{\left (\phi \right )} + 1} - \frac{1}{2 e} \log{\left (\frac{- e \sin{\left (\phi \right )} + 1}{e \sin{\left (\phi \right )} + 1} \right )}\right) \left(\frac{\cos^{2}{\left (\phi_{1} \right )}}{\left(- e^{2} \sin^{2}{\left (\phi_{1} \right )} + 1\right)^{1.0}} - \frac{\cos^{2}{\left (\phi_{2} \right )}}{\left(- e^{2} \sin^{2}{\left (\phi_{2} \right )} + 1\right)^{1.0}}\right)}{- \left(- e^{2} + 1\right) \left(\frac{\sin{\left (\phi_{1} \right )}}{- e^{2} \sin^{2}{\left (\phi_{1} \right )} + 1} - \frac{1}{2 e} \log{\left (\frac{- e \sin{\left (\phi_{1} \right )} + 1}{e \sin{\left (\phi_{1} \right )} + 1} \right )}\right) + \left(- e^{2} + 1\right) \left(\frac{\sin{\left (\phi_{2} \right )}}{- e^{2} \sin^{2}{\left (\phi_{2} \right )} + 1} - \frac{1}{2 e} \log{\left (\frac{- e \sin{\left (\phi_{2} \right )} + 1}{e \sin{\left (\phi_{2} \right )} + 1} \right )}\right)} + \frac{\left(- e^{2} + 1\right) \left(\frac{\sin{\left (\phi_{1} \right )}}{- e^{2} \sin^{2}{\left (\phi_{1} \right )} + 1} - \frac{1}{2 e} \log{\left (\frac{- e \sin{\left (\phi_{1} \right )} + 1}{e \sin{\left (\phi_{1} \right )} + 1} \right )}\right) \left(\frac{\cos^{2}{\left (\phi_{1} \right )}}{\left(- e^{2} \sin^{2}{\left (\phi_{1} \right )} + 1\right)^{1.0}} - \frac{\cos^{2}{\left (\phi_{2} \right )}}{\left(- e^{2} \sin^{2}{\left (\phi_{2} \right )} + 1\right)^{1.0}}\right)}{- \left(- e^{2} + 1\right) \left(\frac{\sin{\left (\phi_{1} \right )}}{- e^{2} \sin^{2}{\left (\phi_{1} \right )} + 1} - \frac{1}{2 e} \log{\left (\frac{- e \sin{\left (\phi_{1} \right )} + 1}{e \sin{\left (\phi_{1} \right )} + 1} \right )}\right) + \left(- e^{2} + 1\right) \left(\frac{\sin{\left (\phi_{2} \right )}}{- e^{2} \sin^{2}{\left (\phi_{2} \right )} + 1} - \frac{1}{2 e} \log{\left (\frac{- e \sin{\left (\phi_{2} \right )} + 1}{e \sin{\left (\phi_{2} \right )} + 1} \right )}\right)} + \frac{\cos^{2}{\left (\phi_{1} \right )}}{\left(- e^{2} \sin^{2}{\left (\phi_{1} \right )} + 1\right)^{1.0}}\right)^{0.5}$$
















In [19]:



    


xe = rhoe * sin(thetae)
ye = rho0e - rhoe * cos(thetae)

I'm not even going to bother printing these two. They're long and it's not important.

Cartopy Boundaries





In [20]:



    


# Default Globe in Cartopy
WGS84_semimajor = Rational(6378137)
WGS84_flattening = Rational(1, 298.257223563)
# Plotting fails if I leave this bit as exact, thinking there are complex values,
# even though the input to sqrt is positive, and it's supposed to return the positive root.
WGS84_ellipticity = sqrt(2 * WGS84_flattening - WGS84_flattening**2).evalf()

# Default parameters in Cartopy
lat0_default = radians(0)
lon0_default = radians(0)
stp1_default = radians(20)
stp2_default = radians(50)

Ellipsoid

Default - WGS84





In [21]:



    


ellipsoid_vars = {
    semimajor: WGS84_semimajor,
    ellipticity: WGS84_ellipticity,
    standard_parallel_1: stp1_default,
    standard_parallel_2: stp2_default,
    central_lat: lat0_default,
    central_lon: lon0_default
}

I know the outer boundary should be something like the following based on existing plots.





In [22]:



    


lat_min_ellipsoid = radians(lat0_default - 90)
lat_max_ellipsoid = radians(lat0_default + 90)
lon_min_ellipsoid = radians(lon0_default - 180)
lon_max_ellipsoid = radians(lon0_default + 180)

default_xe = subdict(xe, ellipsoid_vars)
eq1x = default_xe.subs(point_lat, lat_max_ellipsoid)
eq2x = default_xe.subs(point_lon, lon_max_ellipsoid)
eq3x = default_xe.subs(point_lat, lat_min_ellipsoid)
eq4x = default_xe.subs(point_lon, lon_min_ellipsoid)

default_ye = subdict(ye, ellipsoid_vars)
eq1y = default_ye.subs(point_lat, lat_max_ellipsoid)
eq2y = default_ye.subs(point_lon, lon_max_ellipsoid)
eq3y = default_ye.subs(point_lat, lat_min_ellipsoid)
eq4y = default_ye.subs(point_lon, lon_min_ellipsoid)

plotting.plot_parametric((eq1x, eq1y, (point_lon, lon_min_ellipsoid, lon_max_ellipsoid)),
                         (eq2x, eq2y, (point_lat, lat_max_ellipsoid, lat_min_ellipsoid)),
                         (eq3x, eq3y, (point_lon, lon_max_ellipsoid, lon_min_ellipsoid)),
                         (eq4x, eq4y, (point_lat, lat_min_ellipsoid, lat_max_ellipsoid)),
                         title='Ellipsoid - WGS84 - Guessed Boundaries')



    


















    

















    
Out[22]:








<sympy.plotting.plot.Plot at 0x7f8d25f2e8d0>

















In [23]:



    


lat_magnitude_ellipsoid = subdict((point_lat - central_lat) / radians(90), ellipsoid_vars)
lon_magnitude_ellipsoid = subdict((point_lon - central_lon) / radians(180), ellipsoid_vars)

plotting.plot3d_parametric_surface(default_xe, default_ye, lat_magnitude_ellipsoid,
                                   (point_lat, lat_min_ellipsoid, lat_max_ellipsoid),
                                   (point_lon, lon_min_ellipsoid, lon_max_ellipsoid),
                                   xlabel='X', ylabel='Y',
                                   title='Ellipsoid - WGS84 - Latitude Dependence')
plotting.plot3d_parametric_surface(default_xe, default_ye, lon_magnitude_ellipsoid,
                                   (point_lat, lat_min_ellipsoid, lat_max_ellipsoid),
                                   (point_lon, lon_min_ellipsoid, lon_max_ellipsoid),
                                   xlabel='X', ylabel='Y',
                                   title='Ellipsoid - WGS84 - Longitude Dependence')



    


















    

















    

















    
Out[23]:








<sympy.plotting.plot.Plot at 0x7f8d25ecbda0>

From the above plot, it seems pretty clear that the boundaries match what I guessed above.

However, since the extrema in the $x$ direction is not at the longitudinal limit, we need to take a derivative and find the zeros.





In [24]:



    


xe_deriv = diff(xe, point_lon)
xe_deriv



    


















    
Out[24]:







$$a \left(- \frac{\left(- e^{2} + 1\right) \left(\frac{\sin{\left (\phi \right )}}{- e^{2} \sin^{2}{\left (\phi \right )} + 1} - \frac{1}{2 e} \log{\left (\frac{- e \sin{\left (\phi \right )} + 1}{e \sin{\left (\phi \right )} + 1} \right )}\right) \left(\frac{\cos^{2}{\left (\phi_{1} \right )}}{\left(- e^{2} \sin^{2}{\left (\phi_{1} \right )} + 1\right)^{1.0}} - \frac{\cos^{2}{\left (\phi_{2} \right )}}{\left(- e^{2} \sin^{2}{\left (\phi_{2} \right )} + 1\right)^{1.0}}\right)}{- \left(- e^{2} + 1\right) \left(\frac{\sin{\left (\phi_{1} \right )}}{- e^{2} \sin^{2}{\left (\phi_{1} \right )} + 1} - \frac{1}{2 e} \log{\left (\frac{- e \sin{\left (\phi_{1} \right )} + 1}{e \sin{\left (\phi_{1} \right )} + 1} \right )}\right) + \left(- e^{2} + 1\right) \left(\frac{\sin{\left (\phi_{2} \right )}}{- e^{2} \sin^{2}{\left (\phi_{2} \right )} + 1} - \frac{1}{2 e} \log{\left (\frac{- e \sin{\left (\phi_{2} \right )} + 1}{e \sin{\left (\phi_{2} \right )} + 1} \right )}\right)} + \frac{\left(- e^{2} + 1\right) \left(\frac{\sin{\left (\phi_{1} \right )}}{- e^{2} \sin^{2}{\left (\phi_{1} \right )} + 1} - \frac{1}{2 e} \log{\left (\frac{- e \sin{\left (\phi_{1} \right )} + 1}{e \sin{\left (\phi_{1} \right )} + 1} \right )}\right) \left(\frac{\cos^{2}{\left (\phi_{1} \right )}}{\left(- e^{2} \sin^{2}{\left (\phi_{1} \right )} + 1\right)^{1.0}} - \frac{\cos^{2}{\left (\phi_{2} \right )}}{\left(- e^{2} \sin^{2}{\left (\phi_{2} \right )} + 1\right)^{1.0}}\right)}{- \left(- e^{2} + 1\right) \left(\frac{\sin{\left (\phi_{1} \right )}}{- e^{2} \sin^{2}{\left (\phi_{1} \right )} + 1} - \frac{1}{2 e} \log{\left (\frac{- e \sin{\left (\phi_{1} \right )} + 1}{e \sin{\left (\phi_{1} \right )} + 1} \right )}\right) + \left(- e^{2} + 1\right) \left(\frac{\sin{\left (\phi_{2} \right )}}{- e^{2} \sin^{2}{\left (\phi_{2} \right )} + 1} - \frac{1}{2 e} \log{\left (\frac{- e \sin{\left (\phi_{2} \right )} + 1}{e \sin{\left (\phi_{2} \right )} + 1} \right )}\right)} + \frac{\cos^{2}{\left (\phi_{1} \right )}}{\left(- e^{2} \sin^{2}{\left (\phi_{1} \right )} + 1\right)^{1.0}}\right)^{0.5} \cos{\left (\frac{\left(\lambda - \lambda_{0}\right) \left(\frac{\cos^{2}{\left (\phi_{1} \right )}}{\left(- e^{2} \sin^{2}{\left (\phi_{1} \right )} + 1\right)^{1.0}} - \frac{\cos^{2}{\left (\phi_{2} \right )}}{\left(- e^{2} \sin^{2}{\left (\phi_{2} \right )} + 1\right)^{1.0}}\right)}{- \left(- e^{2} + 1\right) \left(\frac{\sin{\left (\phi_{1} \right )}}{- e^{2} \sin^{2}{\left (\phi_{1} \right )} + 1} - \frac{1}{2 e} \log{\left (\frac{- e \sin{\left (\phi_{1} \right )} + 1}{e \sin{\left (\phi_{1} \right )} + 1} \right )}\right) + \left(- e^{2} + 1\right) \left(\frac{\sin{\left (\phi_{2} \right )}}{- e^{2} \sin^{2}{\left (\phi_{2} \right )} + 1} - \frac{1}{2 e} \log{\left (\frac{- e \sin{\left (\phi_{2} \right )} + 1}{e \sin{\left (\phi_{2} \right )} + 1} \right )}\right)} \right )}$$
















In [25]:



    


xe_ds = subdict(xe_deriv, ellipsoid_vars)
# Known latitude for outer ring
xe_ds = xe_ds.subs(point_lat, radians(-90))



    











In [26]:



    


extrema_ellipsoid = solve(xe_ds, point_lon)
extrema_ellipsoid



    


















    
Out[26]:







$$\left [ 2.83406737943318, \quad 8.50220213829953\right ]$$
















In [27]:



    


extrema_ellipsoid = [degrees(_) for _ in extrema_ellipsoid]
extrema_ellipsoid



    


















    
Out[27]:







$$\left [ 162.380099697222, \quad 487.140299091667\right ]$$
















In [28]:



    


xmin = subdict(xe, ellipsoid_vars)
# Calculated minimum location
xmin = xmin.subs(point_lat, radians(-90)).subs(point_lon, radians(-extrema_ellipsoid[0]))
xmin.evalf()



    


















    
Out[28]:







$$-17702759.799178$$
















In [29]:



    


xmax = subdict(xe, ellipsoid_vars)
# Calculated maximum location
xmax = xmax.subs(point_lat, radians(-90)).subs(point_lon, radians(extrema_ellipsoid[0]))
xmax.evalf()



    


















    
Out[29]:







$$17702759.799178$$
















In [30]:



    


ymin = subdict(ye, ellipsoid_vars)
# Known minimum location from plot
ymin = ymin.subs(point_lat, radians(-90)).subs(point_lon, radians(0))
ymin.evalf()



    


















    
Out[30]:







$$-4782937.05107294$$
















In [31]:



    


ymax = subdict(ye, ellipsoid_vars)
# Known maximum location from plot
ymax = ymax.subs(point_lat, radians(-90)).subs(point_lon, radians(180))
ymax.evalf()



    


















    
Out[31]:







$$15922623.9317694$$

I want to find how many points to split the longitude range into so that it comes really close to this extrema.





In [32]:



    


f = plt.figure(figsize=(16, 16))
for n in range(1, 181, 2):
    lons = np.linspace(-180, 180, n)
    plt.plot(lons, n * np.ones(n), c='b')
    plt.scatter(lons, n * np.ones(n), c='b')
    index = np.argmin(abs(lons - extrema_ellipsoid[0]))
    plt.scatter(lons[index], n, c='r')
plt.axvline(extrema_ellipsoid[0], c='g', ls='--')
plt.xlim(150, 180)
plt.xlabel('Longitude')
plt.ylim(0, 180)
plt.ylabel('Number of Points')
plt.title('Longitude Sampling for a Stipulated Number of Points')
plt.show()

It seems like some good choices are 83, 103, 123, or 165.





In [33]:



    


lons = np.linspace(-180, 180, 103)
index = np.argmin(abs(lons - extrema_ellipsoid[0]))
lons[index-1:index+2], abs(lons[index-1:index+2] - extrema_ellipsoid[0])



    


















    
Out[33]:








(array([ 158.82352941,  162.35294118,  165.88235294]),
 array([3.55657028545752, 0.0271585207516694, 3.50225324395424], dtype=object))

Eccentric Globe





In [34]:



    


a = Rational(1000)
b = Rational(500)
f = (a - b) / a
eccentric_globe_vars = {
    semimajor: a,
    ellipticity: sqrt(2 * f - f**2).evalf(),
    standard_parallel_1: stp1_default,
    standard_parallel_2: stp2_default,
    central_lat: lat0_default,
    central_lon: lon0_default
}



    











In [35]:



    


lat_min_eccentric_globe = radians(lat0_default - 90)
lat_max_eccentric_globe = radians(lat0_default + 90)
lon_min_eccentric_globe = radians(lon0_default - 180)
lon_max_eccentric_globe = radians(lon0_default + 180)

eccentric_xe = subdict(xe, eccentric_globe_vars)
eq1x = eccentric_xe.subs(point_lat, lat_max_eccentric_globe)
eq2x = eccentric_xe.subs(point_lon, lon_max_eccentric_globe)
eq3x = eccentric_xe.subs(point_lat, lat_min_eccentric_globe)
eq4x = eccentric_xe.subs(point_lon, lon_min_eccentric_globe)

eccentric_ye = subdict(ye, eccentric_globe_vars)
eq1y = eccentric_ye.subs(point_lat, lat_max_eccentric_globe)
eq2y = eccentric_ye.subs(point_lon, lon_max_eccentric_globe)
eq3y = eccentric_ye.subs(point_lat, lat_min_eccentric_globe)
eq4y = eccentric_ye.subs(point_lon, lon_min_eccentric_globe)

plotting.plot_parametric((eq1x, eq1y, (point_lon, lon_min_eccentric_globe, lon_max_eccentric_globe)),
                         (eq2x, eq2y, (point_lat, lat_max_eccentric_globe, lat_min_eccentric_globe)),
                         (eq3x, eq3y, (point_lon, lon_max_eccentric_globe, lon_min_eccentric_globe)),
                         (eq4x, eq4y, (point_lat, lat_min_eccentric_globe, lat_max_eccentric_globe)),
                         title='Ellipsoid - Eccentric Globe - Guessed Boundaries')



    


















    

















    
Out[35]:








<sympy.plotting.plot.Plot at 0x7f8d1fa85eb8>

















In [36]:



    


lat_magnitude_eccentric_globe = subdict((point_lat - central_lat) / radians(90), eccentric_globe_vars)
lon_magnitude_eccentric_globe = subdict((point_lon - central_lon) / radians(180), eccentric_globe_vars)

plotting.plot3d_parametric_surface(eccentric_xe, eccentric_ye, lat_magnitude_eccentric_globe,
                                   (point_lat, lat_min_eccentric_globe, lat_max_eccentric_globe),
                                   (point_lon, lon_min_eccentric_globe, lon_max_eccentric_globe),
                                   xlabel='X', ylabel='Y',
                                   title='Ellipsoid - Eccentric Globe - Latitude Dependence')
plotting.plot3d_parametric_surface(eccentric_xe, eccentric_ye, lon_magnitude_ellipsoid,
                                   (point_lat, lat_min_eccentric_globe, lat_max_eccentric_globe),
                                   (point_lon, lon_min_eccentric_globe, lon_max_eccentric_globe),
                                   xlabel='X', ylabel='Y',
                                   title='Ellipsoid - Eccentric Globe - Longitude Dependence')



    


















    

















    

















    
Out[36]:








<sympy.plotting.plot.Plot at 0x7f8d2b0ce780>

It appears for the eccentric globe, a very similar boundary will occur.





In [37]:



    


xe_deg = subdict(xe_deriv, ellipsoid_vars)
# Known latitude for outer ring
xe_deg = xe_deg.subs(point_lat, radians(-90))



    











In [38]:



    


extrema_eccentric_globe = solve(xe_deg, point_lon)
extrema_eccentric_globe



    


















    
Out[38]:







$$\left [ 2.83406737943318, \quad 8.50220213829953\right ]$$
















In [39]:



    


extrema_eccentric_globe = [degrees(_) for _ in extrema_eccentric_globe]
extrema_eccentric_globe



    


















    
Out[39]:







$$\left [ 162.380099697222, \quad 487.140299091667\right ]$$

These are (unsurprisingly) at the same place!





In [40]:



    


xmin = subdict(xe, eccentric_globe_vars)
# Calculated minimum location
xmin = xmin.subs(point_lat, radians(-90)).subs(point_lon, radians(-extrema_eccentric_globe[0]))
xmin.evalf()



    


















    
Out[40]:







$$-2323.47073363411$$
















In [41]:



    


xmax = subdict(xe, eccentric_globe_vars)
# Calculated maximum location
xmax = xmax.subs(point_lat, radians(-90)).subs(point_lon, radians(extrema_eccentric_globe[0]))
xmax.evalf()



    


















    
Out[41]:







$$2323.47073363411$$
















In [42]:



    


ymin = subdict(ye, eccentric_globe_vars)
# Known minimum location from plot
ymin = ymin.subs(point_lat, radians(-90)).subs(point_lon, radians(0))
ymin.evalf()



    


















    
Out[42]:







$$-572.556243423971$$
















In [43]:



    


ymax = subdict(ye, eccentric_globe_vars)
# Known maximum location from plot
ymax = ymax.subs(point_lat, radians(-90)).subs(point_lon, radians(180))
ymax.evalf()



    


















    
Out[43]:







$$2402.36176984391$$
















In [44]:



    


ymax = subdict(ye, eccentric_globe_vars)
# Known maximum location from plot
ymax = ymax.subs(point_lat, radians(-90)).subs(point_lon, radians(180))
ymax.evalf()



    


















    
Out[44]:







$$2402.36176984391$$

Clarke 1866





In [45]:



    


# See example on pg 292
clarke1866_lat0 = 23
clarke1866_lon0 = -96
clarke1866_vars = {
    semimajor: 6378206.4,
    ellipticity: 0.0822719,
    standard_parallel_1: radians(29 + 30 / 60),
    standard_parallel_2: radians(45 + 30 / 60),
    central_lat: radians(clarke1866_lat0),
    central_lon: radians(clarke1866_lon0)
}



    











In [46]:



    


lat_min_clarke1866 = radians(-90)
lat_max_clarke1866 = radians(+90)
lon_min_clarke1866 = radians(clarke1866_lon0 - 180)
lon_max_clarke1866 = radians(clarke1866_lon0 + 180)

clarke1866_xe = subdict(xe, clarke1866_vars)
eq1x = clarke1866_xe.subs(point_lat, lat_max_clarke1866)
eq2x = clarke1866_xe.subs(point_lon, lon_max_clarke1866)
eq3x = clarke1866_xe.subs(point_lat, lat_min_clarke1866)
eq4x = clarke1866_xe.subs(point_lon, lon_min_clarke1866)

clarke1866_ye = subdict(ye, clarke1866_vars)
eq1y = clarke1866_ye.subs(point_lat, lat_max_clarke1866)
eq2y = clarke1866_ye.subs(point_lon, lon_max_clarke1866)
eq3y = clarke1866_ye.subs(point_lat, lat_min_clarke1866)
eq4y = clarke1866_ye.subs(point_lon, lon_min_clarke1866)

plotting.plot_parametric((eq1x, eq1y, (point_lon, lon_min_clarke1866, lon_max_clarke1866)),
                         (eq2x, eq2y, (point_lat, lat_max_clarke1866, lat_min_clarke1866)),
                         (eq3x, eq3y, (point_lon, lon_max_clarke1866, lon_min_clarke1866)),
                         (eq4x, eq4y, (point_lat, lat_min_clarke1866, lat_max_clarke1866)),
                         title='Ellipsoid - Clarke 1866 - Guessed Boundaries')



    


















    

















    
Out[46]:








<sympy.plotting.plot.Plot at 0x7f8d25b08be0>

















In [47]:



    


lat_magnitude_clarke1866 = subdict((point_lat - central_lat) / radians(90), clarke1866_vars)
lon_magnitude_clarke1866 = subdict((point_lon - central_lon) / radians(180), clarke1866_vars)

plotting.plot3d_parametric_surface(clarke1866_xe, clarke1866_ye, lat_magnitude_clarke1866,
                                   (point_lat, lat_min_clarke1866, lat_max_clarke1866),
                                   (point_lon, lon_min_clarke1866, lon_max_clarke1866),
                                   xlabel='X', ylabel='Y',
                                   title='Ellipsoid - Clarke 1866 - Latitude Dependence')
plotting.plot3d_parametric_surface(clarke1866_xe, clarke1866_ye, lon_magnitude_clarke1866,
                                   (point_lat, lat_min_clarke1866, lat_max_clarke1866),
                                   (point_lon, lon_min_clarke1866, lon_max_clarke1866),
                                   xlabel='X', ylabel='Y',
                                   title='Ellipsoid - Clarke 1866 - Longitude Dependence')



    


















    

















    

















    
Out[47]:








<sympy.plotting.plot.Plot at 0x7f8d25aeea58>

It appears for Clarke 1866, a very similar boundary will occur once again.





In [48]:



    


xe_dc = subdict(xe_deriv, clarke1866_vars)
# Known latitude for outer ring
xe_dc = xe_dc.subs(point_lat, radians(-90))



    











In [49]:



    


extrema_clarke1866 = solve(xe_dc, point_lon)
extrema_clarke1866



    


















    
Out[49]:







$$\left [ 0.92986989604966, \quad 6.14064185197808\right ]$$
















In [50]:



    


extrema_clarke1866 = [degrees(_) for _ in extrema_clarke1866]
extrema_clarke1866



    


















    
Out[50]:







$$\left [ 53.2776205399141, \quad 351.832861619742\right ]$$
















In [51]:



    


xmin = subdict(xe, clarke1866_vars)
# Calculated minimum location
xmin = xmin.subs(point_lat, radians(-90)).subs(point_lon, radians(extrema_clarke1866[1]))
xmin.evalf()



    


















    
Out[51]:







$$-16900972.674607$$
















In [52]:



    


xmax = subdict(xe, clarke1866_vars)
# Calculated maximum location
xmax = xmax.subs(point_lat, radians(-90)).subs(point_lon, radians(extrema_clarke1866[0]))
xmax.evalf()



    


















    
Out[52]:







$$16900972.674607$$
















In [53]:



    


ymin = subdict(ye, clarke1866_vars)
# Known minimum location from plot
ymin = ymin.subs(point_lat, radians(-90)).subs(point_lon, radians(clarke1866_lon0))
ymin.evalf()



    


















    
Out[53]:







$$-6971893.11311231$$
















In [54]:



    


ymax = subdict(ye, clarke1866_vars)
# Known maximum location from plot
ymax = ymax.subs(point_lat, radians(-90)).subs(point_lon, lon_max_clarke1866)
ymax.evalf()



    


















    
Out[54]:







$$15298166.8919989$$

Sphere





In [55]:



    


# From numerical examples, pg 291
radius_sphere = 1.0
stp1_sphere = 29 + 30 / 60
stp2_sphere = 45 + 30 / 60
lat0_sphere = 23
lon0_sphere = -96

sphere_vars = {radius: radius_sphere,
               standard_parallel_1: radians(stp1_sphere),
               standard_parallel_2: radians(stp2_sphere),
               central_lat: radians(lat0_sphere),
               central_lon: radians(lon0_sphere)}



    











In [56]:



    


lon_min_sphere = radians(lon0_sphere - 180)
lon_max_sphere = radians(lon0_sphere + 180)
# Interestingly, the boundary does not appear to depend on central latitude!
lat_min_sphere = radians(-90)
lat_max_sphere = radians(+90)

sphere_xs = subdict(xs, sphere_vars)
eq1x = sphere_xs.subs(point_lat, lat_max_sphere)
eq2x = sphere_xs.subs(point_lon, lon_max_sphere)
eq3x = sphere_xs.subs(point_lat, lat_min_sphere)
eq4x = sphere_xs.subs(point_lon, lon_min_sphere)

sphere_ys = subdict(ys, sphere_vars)
eq1y = sphere_ys.subs(point_lat, lat_max_sphere)
eq2y = sphere_ys.subs(point_lon, lon_max_sphere)
eq3y = sphere_ys.subs(point_lat, lat_min_sphere)
eq4y = sphere_ys.subs(point_lon, lon_min_sphere)

plotting.plot_parametric((eq1x, eq1y, (point_lon, lon_min_sphere, lon_max_sphere)),
                         (eq2x, eq2y, (point_lat, lat_max_sphere, lat_min_sphere)),
                         (eq3x, eq3y, (point_lon, lon_max_sphere, lon_min_sphere)),
                         (eq4x, eq4y, (point_lat, lat_min_sphere, lat_max_sphere)),
                         title='Sphere - Guessed Boundaries')



    


















    

















    
Out[56]:








<sympy.plotting.plot.Plot at 0x7f8d1ffb5e80>

Evidently, the boundary does not depend on central latitude. However, central longitude does seem to make a difference.





In [57]:



    


lat_magnitude = subdict((point_lat - central_lat) / radians(90), sphere_vars)
lon_magnitude = subdict((point_lon - central_lon) / radians(180), sphere_vars)

plotting.plot3d_parametric_surface(sphere_xs, sphere_ys, lat_magnitude,
                                   (point_lat, lat_min_sphere, lat_max_sphere),
                                   (point_lon, lon_min_sphere, lon_max_sphere),
                                   xlabel='X', ylabel='Y',
                                   title='Sphere - Latitude Dependence')

plotting.plot3d_parametric_surface(sphere_xs, sphere_ys, lon_magnitude,
                                   (point_lat, lat_min_sphere, lat_max_sphere),
                                   (point_lon, lon_min_sphere, lon_max_sphere),
                                   xlabel='X', ylabel='Y',
                                   title='Sphere - Longitude Dependence')



    


















    

















    

















    
Out[57]:








<sympy.plotting.plot.Plot at 0x7f8d240537b8>

The first plot above points to the location of the minimum and maximum latitudes. The second plot is for the longitudes. Clearly the outer boundary occurs at the minimum of the latitude ($-90^\circ$!), and all along the longitudes.

Since the extrema in the $x$ direction is not at the longitudinal limit, we need to take a derivative and find the zeros.





In [58]:



    


xs_deriv = diff(xs, point_lon)
xs_deriv



    


















    
Out[58]:







$$R \left(- \left(\sin{\left (\phi_{1} \right )} + \sin{\left (\phi_{2} \right )}\right) \sin{\left (\phi \right )} + \left(\sin{\left (\phi_{1} \right )} + \sin{\left (\phi_{2} \right )}\right) \sin{\left (\phi_{1} \right )} + \cos^{2}{\left (\phi_{1} \right )}\right)^{0.5} \cos{\left (\left(\lambda - \lambda_{0}\right) \left(\frac{1}{2} \sin{\left (\phi_{1} \right )} + \frac{1}{2} \sin{\left (\phi_{2} \right )}\right) \right )}$$
















In [59]:



    


xs_ds = subdict(xs_deriv, sphere_vars)
# Known latitude for outer ring
xs_ds = xs_ds.subs(point_lat, lat_min_sphere)
xs_ds



    


















    
Out[59]:







$$1.59902949775029 \cos{\left (0.602837004628824 \lambda + 1.0100630960288 \right )}$$
















In [60]:



    


extrema_sphere = solve(xs_ds, point_lon)
extrema_sphere



    


















    
Out[60]:







$$\left [ 0.9301572837443, \quad 6.14150401506203\right ]$$
















In [61]:



    


extrema_sphere = [degrees(_) for _ in extrema_sphere]
extrema_sphere



    


















    
Out[61]:







$$\left [ 53.2940866419009, \quad 351.882259925704\right ]$$
















In [62]:



    


xmin = subdict(xs, sphere_vars)
# Calculated minimum location
xmin = xmin.subs(point_lat, lat_min_sphere).subs(point_lon, radians(-extrema_sphere[0]))
xmin.evalf()



    


















    
Out[62]:







$$1.15214837563775$$
















In [63]:



    


xmax = subdict(xs, sphere_vars)
# Calculated maximum location
xmax = xmax.subs(point_lat, lat_min_sphere).subs(point_lon, radians(extrema_sphere[0]))
xmax.evalf()



    


















    
Out[63]:







$$2.6525072042232$$

Obviously, something has gone wrong here. Maybe by trying the other extrema, the result will make some sense.





In [64]:



    


xmin = subdict(xs, sphere_vars)
# Calculated minimum location
xmin = xmin.subs(point_lat, lat_min_sphere).subs(point_lon, radians(-extrema_sphere[1]))
xmin.evalf()



    


















    
Out[64]:







$$-1.15214837563774$$
















In [65]:



    


xmax = subdict(xs, sphere_vars)
# Calculated maximum location
xmax = xmax.subs(point_lat, lat_min_sphere).subs(point_lon, radians(extrema_sphere[1]))
xmax.evalf()



    


















    
Out[65]:







$$-2.6525072042232$$

Unlike the ellipsoid cases above, it appears both roots were needed to see the full limits. This difference is likely due to the central longitude not being zero, and therefore the roots are not symmetric.





In [66]:



    


ymin = subdict(ys, sphere_vars)
# Known minimum location from plot
ymin = ymin.subs(point_lat, lat_min_sphere).subs(point_lon, radians(lon0_sphere))
ymin.evalf()



    


















    
Out[66]:







$$-1.09628087472359$$
















In [67]:



    


ymax = subdict(ys, sphere_vars)
# Known maximum location from plot
ymax = ymax.subs(point_lat, lat_min_sphere).subs(point_lon, lon_max_sphere)
ymax.evalf()



    


















    
Out[67]:







$$2.3983472405755$$

Content source: QuLogic/map-projections

Similar notebooks:

notebook.community | gallery | about